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Abstract 

Hybrid QM(DFT) /MM molecular dynamics simulations have been carried out for the Watson- 
Crick base pair of 9-ethyl-8-phenyladenine and 1-cyclohexyluracil in deutero chloroform solution 
at room temperature. Trajectories are analyzed putting special attention to the geometric 
correlations of the N-H • • • N and N-H • • • O hydrogen bonds in the base pair. Further, based 
on empirical correlations between the hydrogen bond bond length and the fundamental NH 
stretching frequency its fluctuations are obtained along the trajectory. Using the time dependent 
frequencies the infrared lineshape is determined assuming the validity of a second order cumulant 
expansion. The deviations for the fundamental transition frequencies are calculated to amount 
to less than 2% as compared with experiment. The width of the spectrum for the N-H • • • N bond 
is in reasonable agreement with experiment while that for the N-H • • • O case is underestimated 
by the present model. Comparing the performance of different pseudopotentials it is found that 
the Troullier-Martins pseudopotential with a 70 Ry cut-off yields the best agreement. 

PACS numbers: 71.15.Pd, 82.30.Rs, 82.39.Pj, 87.16.dt, 87.64.km 
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I. INTRODUCTION 



One of the primary examples for the importance of hydrogen bonding is its role in the 
structural selectivity of the pairing of nucleic acid bases in DNA Here double and 
triple hydrogen bonds (HBs) are formed whose theoretical characterization still poses a 
challenge as it combines many features of multidimensional condensed phase quantum 
dynamics. Infrared (IR) spectra are particularly sensitive to HB formation. The broad 
lineshapes not only contain information on the intramolecular couplings, but the underly- 
ing fluctuation dynamics also reflects the interaction with the surrounding medium [2], yl . 
The spectroscopic signatures of hydrogen bonding in heterodimers of various bases have 
been investigated in different phases. Isolated adenine:thymine pairs, for instance, were 
studied in gas phase using IR-UV double resonance experiments j4|. This provided evi- 
dence that Watson-Crick pairing does not result in the most probable form under these 
conditions. Later, the dominant tautomer could be identifled on the basis of quantum 
dynamical simulations of IR spectra including anharmonicity Ultrafast nonlinear 
IR spectroscopy has been used to scrutinize the vibrational dynamics of individual base 
pairs in solution j^. Here it was found, for instance, that the life time of the H-bonded 
NH stretching fundamental transition decreased by a factor of three as compared to the 
monomer case, thus demonstrating the effect of HB mediated anharmonic couplings. The 
next step in complexity has been addressed by ultrafast two-color IR studies of the HBs 
in DNA oligomers. Heyne et al. showed that excitation in the flngerprint region provides, 
by virtue of the anharmonic couplings, a means to identify, for instance, the symmetric 
NH2 stretching vibration of adenine:thymine pairs upon probing in the 3200 cm~^ range 







In the linear IR spectrum this range is completely masked by the absorption of water 
molecules. Subsequently, the assignment has been conflrmed by anisotropy measurements 

a. 

Since flrst principles quantum dynamical simulation of condensed phase HB dynamics 
will remain elusive different approximate strategies are usually followed. In order to cap- 
ture the electronic structure associated with H-bonding density functional theory (DFT) 
seems to be most appropriate as long as stacking interactions are of minor importance 
jol. By running classical trajectories with on-the-fly DFT forces the IR spectrum can be 
obtained from the Fourier transform of the dipole autocorrelation function. Of course, 
this completely neglects quantum effects in the nuclear motion. The latter are frequently 
accounted for by quantum correction factors [l^. Car-Parrinello molecular dynamics 



(CPMD) has proven to be particularly well-suited for simulating condensed phase dynam- 
ics including H-bonded systems. A pioneering contribution has been given by Silvestrelli 
et al. who determined the IR spectrum of liquid water [ll|. Subsequent applications to 
aqueous solutions, e.g., of uracil have been summarized in Ref. jl^]- More recent studies 
focussed on strongly H-bonded system like the Zundel cation in the crystalline phase jisl . 

In principle quantum effects for the fast proton stretching motion can be recovered by 
calculating its potential energy curve for selected configurations along the CPMD trajec- 
tory (snapshot potentials). Upon solving the stationary Schrodinger equation for these 
potentials one obtains a set of fundamental transition frequencies. The combination of 
CPMD and snapshot potentials for the proton stretching motion has been put forward in 
Ref. jl^ . In this reference an intramolecular N-H ■ ■ ■ O HB in a Mannich base in gas and 
crystalline phase was considered. The IR spectrum has been calculated by convolution 
of the stick spectrum obtained along the trajectory with a Gaussian. This gave a band 
maximum below 2000 cm~^ , considerably red shifted with respect to the estimated exper- 
imental value. Interestingly, the analogous simulation in CCI4 solution gave a band whose 
first moment has been well above 2000 cm~^ , pointing out the sensitivity with respect to 
the environment In Ref. 16 1 the method was applied to an A^-oxide which exhibits 
two rather strong nonequivalent 0-H ■ ■ ■ O HBs in the unit cell. Here, two-dimensional 
snapshot potentials including the bending vibration of the HB had been considered as 
well. Finally, a different A^-oxide in the crystalline phase with a rather strong HB was 
studied in Ref. In this case a 1000 cm~^ broad distribution of transition frequencies 
for the OH-stretching fundamental was found from one- dimensional snapshot potentials 
in reasonable agreement with experiment. 

Quantum Mechanics/Molecular Mechanics (QM/MM) hybrid methods have been de- 
veloped to cope with situations where only a part of a large system needs to be treated at 
the quantum mechanical level (for a review, see Ref. [isl). Here the crucial point is the 
interface between the QM and MM parts and in the context of CPMD different schemes 
for the most delicate issue of the nonbonding electrostatic interactions have been devel- 
oped [l9|,0. A series of impressive applications to the IR spectroscopy of various forms 



of water clusters in a model of bacteriorhodopsin has been presented in Refs. |2ll . |22| . |23 |. 

Recently, we have applied the QM/MM approach to the determination of the IR line- 
shape of the uracil NH stretching vibration in a modified adenine: uracil (A:U) pair in 
deuterochloroform solution [24 1. In doing so we made use of an empirical NH-frequency- 
N ■ ■ ■ N-distance correlation which gave access to the fluctuating transition frequency along 



the trajectory (see also Ref. 25|). From the autocorrelation function of this transition 
frequency the lineshape could be determined. This procedure gave the band maximum at 
3209 cm~^ and the full-width at half maximum (FWHM) of 39 cm~^ , both values being 

n 

in reasonable agreement with experiment j6|. 

In the present contribution we will extend our previous study in several respects. 
First, the IR spectral features of both intermolecular HBs are investigated. This requires 
to establish a distance-frequency correlation for N-H ■ ■ ■ O HBs which is done on the 
basis of available crystal structure data. Alternatively, we explore the possibility to use 
snapshot potentials for obtaining this correlation curve. It will turn out that the snapshot 
frequencies underestimate the average transition frequency considerably and are of little 
use for the considered intermolecular HBs. In order to scrutinize this issue we report 
an investigation of the influence on the used pseudopotential (PP) and kinetic energy 
cut-off parameter in the DFT calculation. Besides the IR spectrum we will put emphasis 
on the geometric correlations across the double HBs which can be obtained along the 
trajectory. This concerns in particular the linearity, planarity, and the correlation between 
HB distance and proton position. 

The paper is organized as follows: In the next Section II we start by defining the 
model system and the QM/MM protocol. Further, the issues of geometric correlations 
and lineshape analysis are discussed. In Section HI we first give results on geometric 
correlations before we come to the lineshape and its dependence on the PP used. The 
results are summarized in Section IV. 



II. METHODS 

A. QM/MM Protocol 

We will investigate the vibrational dynamics of 9-ethyl-8-phenyladenine (A) and 
1-cyclohexyluracil (U) solvated in CDCI3 which has been studied experimentally by 
Woutersen and Cristalli [6j. To this end a hybrid QM/MM method is used as implemented 
in the Gromacs/CPMD interface [2^. All atoms of the A:U pair except those belonging 
to the substituents are dealt with quantum mechanically by the CPMD package [26]. The 
substituents themselves as well as the solvent are described by the OPLS all-atom force 
field as implemented in Gromacs [28]. H-atom capping is used to saturate the dan- 
gling bonds at the QM/MM boundary (see Fig. [1]). In the simulation a single A:U pair is 




FIG. 1: QM part of the A:U base pair with the capping atoms indicated by circles. In the 
analysis below we will use the angles a (Ng-04-C4) and /5 (C4-N3-N1), the angle 7 between the 
two vectors N3N1 and O4N6 as well the dihedral angle 4> defined by the atoms Ni, N3, O4, and 
Ne. 

solvated in 100 CDCI3 molecules within a box having dimensions 30.0 A x 23.5 A x 23.5 A 
(density is 0.1 M [6]). The QM part is placed in a 21.2 A x 15.9 A x 15.9 A box. The Becke 
exchange and Lee- Yang-Par correlation functional (BLYP) together with the plane wave 
basis set is used as implemented in the CPMD code. The MM molecular dynamics run is 
performed at 298 K using a time step of 2 fs. Note that in the current Gromacs/CPMD 
implementation it is actually a Born-Oppenheimer simulation which is performed. For the 
initial configuration the geometry of the optimized gas phase base pair replacing solvent 
atoms in the equilibrated simulation box has been used. First, a 1 ps equilibration run 
was performed using the Vanderbilt (VB) ultrasoft PP with a plane wave cutoff of 30 Ry. 
Subsequently, trajectories were propagated up to 5.0 ps using different PPs and cut-off 



parameters, i.e., VB 



29] with a cut-off of 30 Ry and 40 Ry, Troullier-Martins (TM) PP 
30( with 70 Ry and 100 Ry, and Goedecker (GO) PP [31] with 100 Ry and 140 Ry. 
Except for the comparative studies the results presented are based on the TM PP using 
a cut-off of 70 Ry. 



B. Hydrogen Bond Geometry 

Hydrogen bonds A-H- ■ • B display remarkable correlations in their geometry. Limbach 
and coworkers promoted a model which gives a simple explanation of these correlation in 
terms of valence bond orders viewing the HB as composed of two diatomic units A-H 



and B-H with bond lengths R\ and i?2, respectively |32|]: 



Pi = exp{-(i?i - Ft°')lhi) , 



(1) 
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FIG. 2: Geometric parameters of the HBs along the trajectory: (a) and (b) HB lengths of 
N-H • • • N and N-H • • • O, respectively. The horizontal line indicates the time average, (c) and 
(d) out-of-line motion of the H atom in the two HBs, measured by the difference L(N-H) + 
L(N • • • H) - L(N • • • N) and L(N-H) + L(0 • • • H) - L(N • • • O), respectively. 



where R^"^ is the equihbrium bond length of the hypothetical nonbonded diatomic and bt 
is a parameter describing the change of the bond valence upon bond stretching. If one 
assumes that the total bond order is unity, i.e. P1+P2 = 1, the two bond lengths cannot 
be changed independently. This can be alternatively expressed in terms of the deviation 
from the HB center — -R2)/2 and the HB length Ri + R2. The correlation curve 
so obtained from numerous structural data expresses a fact well known from quantum 
chemical studies of potential energy surfaces, namely that the HB is compressed upon H 
transfer while passing the transition state [2]. This correlation for a given type of HB - 
although expressed in Eq. ([1]) by four parameters only - yields a rather robust description 



of the HB geometries (see, e.g. the case of N-H- ■ -N bonds in Refs. 32, ISSj). It turns 



out, however, that in particular for short HBs the original formulation required some 



modification. In Ref. 



34| Limbach et al. suggested an empirical correction as follows 



Pl/2 = Pl/2 T C{pi - P2){piP2f - d{piP2f 



(2) 



for which pf + p^ < 1 holds. Here, c and d are parameters to be determined by fitting 
geometric correlations to experimental data. 

The change in the potential energy surface which is reflected in the above geometric 



correlation also manifests itself in a dependence of the AH stretching frequency on the 
HB length. In a seminal work Novak collected experimental data on different H-bonded 
crystals to establish A-B bond length-frequency correlation for N-H ■ ■ ■ N and 0-H- ■ ■ O 



bonds 



351 ] ■ Later, this study was extended to asymmetric HBs by Mikenda [36]. For the 



case of N-H- ■ ■ N bonds the correlation has been expressed previously by the function 24] 



/(r) = cUoo/27rcerf a^j , (3) 

where cUqo is the free NH stretching frequency in gas phase (for parameters see caption 
of Fig. [6]). Eq. ([3]) can also be used for asymmetric HBs such as those of N-H ■ ■ ■ O 
type involving also HBs formed by amine groups (see Supplementary Information). The 
advantage of such relation is obvious, it allows to assign stretching frequencies along a 



trajectory on the basis of atomic coordinates only |24| . 



C. IR Lineshape 

In the classical limit the IR lineshape can be calculated straightforwardly from the 
trajectory by means of the dipole-dipole correlation function 45j. Here we will use line- 
shape theory and account for the quantum nature of the high-frequency proton stretching 
vibrations. Specifically, we will calculate their gap correlation function, defined for some 
fundamental stretching frequency, uJio{t), along a trajectory as follows: 

C{t) = {{u.oit) - {u.o)) (a;io(0) - {uio))) , (4) 

where (cuio) denotes the mean value along the trajectory. Assuming a second order cu- 
mulant expansion to hold the IR spectrum can be expressed via the lineshape function 

g{t)^ f dr rdr'Cir'), (5) 
Jo Jo 



37] 

1 POO 

a{u) = -Re ^^e'("-<"^°>)*-^(*). (6) 
Jo 

Of course, Eq. (jlj), is a classical quantity and does not satisfy detailed balance. There 
are various recipes to compensate for this deficiency by introducing quantum correction 
factors fl^. Recently, Skinner and coworker gave an interesting account on this issue 
which revealed that for the exemplary case of HOD in D2O quantum corrections have a 
only modest influence on the IR lineshape 381] . In order to take into account quantum 



effects we will fit the classical correlation function to the expression obtained within the 



multimode oscillator model 



37| 



C{t) = S{ujj)uj'^ cot\i(hujj /2k-BT) cos{uj jt) 
j 

+ i S{ujj)u^ sm{uj jt) , (7) 
j 

where S{ujj) is the dimensionless Huang- Rhys factor giving the coupling strength between 
the two-level system and the bath mode with frequency ujj. These two sets of parameters 
are used for fitting the real part of C{t) obtained from the classical molecular dynamics 
which in turn gives the imaginary part as well. 

III. RESULTS 

A. Geometry-based Correlations 

In Fig. [2^ and [2]d we depict the N ■ ■ ■ N and the N ■ ■ ■ O distance, respectively, along 
the production part of the trajectory. In the course of the trajectory no H-atom transfer 
occurs. The average HB lengths are L(N ■ ■ ■ N) = 3.09 A and L(N ■ ■ ■ O) = 3.03 A. Since 
geometry-based correlations are usually discussed assuming linear HBs we illustrate in 
panels (c) and (d) of Fig. [2] the deviation of the H atoms from linear motion. The average 
deviations are 0.028 A for the N-H ■ ■ ■ N and 0.038 A for the N-H ■ ■ • O bond, that is, both 
cases can be considered to correspond to almost linear HBs. 

The planarity of the gas phase A:U heterocylces is no longer present in solution. Fo- 
cussing on the two HBs only, one can define the dihedral angle formed by the atoms Ni, 
N3, O4, and Ng (see. Fig. [T]). Its average value during the propagation interval is -19°. 
Fig. [3] shows how this deformation is related to other geometrical changes of the HBs. 
From panel (a) we see that nonplanarity comes along with a decrease of the angle a + f3, 
that is, with a compression of the HBs due to N3HU and/or C4O4 bending vibrations. Fig. 
[3]d reveals an almost perfect linear correlation between the dihedral angle and the angle 
7 describing the directions of the two HBs. Besides the geometric correlation between 
two HBs, there also exists correlation for each HB separately. Fig. H] gives a test of the 



empirical correlation between N-H and HB length derived from 
from the robustness of the correlation reported before 



32 



33 



)q. As anticipated 



34| . geometries along the 



trajectory are remarkably well described by this simple expression. This holds true, irre- 
spective of the deviation from linearity and planarity mentioned above. Moreover, there 
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FIG. 3: Geometric correlations of the N-H • • • N and N-H • • • O HBs along the QM/MM trajectory 
between the absolute value of the dihedral angle (p and (a) the sum of the angles a and /3 and 
(b) the angle 7 (for definitions see caption of Fig. [T]). 
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FIG. 4: Geometric correlations between HB length and H-atom position following from a fit 
to Eq. ([2]) (solid line) with the points sampled along the QM/MM trajectory, (a) N-H • • • N 
(1 = 2 = NH) with Rl"^ = 1.025 A, h = 0.361 A, c = 577.64, and d = 0.30 and (b) N-H • • • O 
(1 = NH, 2 = OH) with iif = 1.021 A, h = 0.407 A, Rl"^ = 1.029 A, 62 = 0.232 A, c = 400, 



d= 0.21. In panel (a) we also show the correlation curve from Limbach et al. (dashed line) 



2m- 



seems to be little variance with respect to the class of molecules as shown by plotting the 
correlation curve with parameters obtained by Limbach et al. for intramolecular N-H ■ ■ ■ N 
HBs (dashed line in panel (a)) 33 1. 

In a next step we explore the correlation between HB bond lengths N ■ ■ ■ N and N ■ ■ ■ O 
and the fundamental transition frequencies of the NH stretching motion. A straightfor- 
ward realization would consist in the calculation of NH stretching snapshot potentials in 




L(N-H)[Ang.] 

FIG. 5: Comparison of potential energy curves for the N-H • • • N HB as calculated from different 
PPs by changing the N-H bond length (11 points in the given interval) with all other coordinates 
being fixed at the t = ps structure. The horizontal lines correspond to the lowest vibrational 
states which have been obtained using the Fourier Grid Hamiltonian method [39]. The cut-off 
parameters are 30 Ry for VB, 70 Ry for TM and 140 Ry for GO. 

an otherwise frozen geometry for some representative points along the trajectory. Inter- 
polating these data one will have at hand a distance-frequency correlation to be used for 
all points along the trajectory. Fig. [5] depicts a representative potential for the N-H ■ ■ ■ N 
HB and its dependence on the form of the PP. First, we notice that up to the overtone 
excitation the potential is rather well described by an anharmonic oscillator and the sec- 
ond rather shallow minimum corresponding to an H transfer configuration is energetically 
not favorable. Second, apart from the shape of this shallow transfer minimum the three 
PPs give rather similar fundamental transition frequencies, which are 2634 cm~^ for VB, 
2647 cm~^ for TM and 2600 for GO. Third, the fundamental transition frequency in the 
2600 cm~^ range is considerable lower than the experimental value of 3185 cm~^ |6|. 

A more complete picture is obtained from Fig. [6^ where we show results from snapshot 
frequency calculations with the VB, TM and GO PPs sampled along the trajectory. They 
reproduce the expected qualitative dependence on the HB length, but are red-shifted by 
about 450 cm~^ with respect to the empirical correlation curve obtained from crystal data. 
This holds irrespective of the PP and cut-off parameter. 

The first question to be clarified concerns the experimental assignment in Ref. j^. 
In fact the similar case of 9-ethyladenine and cyclohexyluracil in chloroform has been 
studied previously and it was found that self-association of U:U dimers occurs which 
would have an NH stretching absorption in the 3000-3200 cm^^ range as well. However, 
the association constant for U:U dimers is 15 times smaller than for A:U ones, so that 
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FIG. 6: (a) Correlation between HB length N---N and fundamental NH transition fre- 
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quency (bullets: spectroscopic data for crystals containing intermolecular N-H • • • N HB 
QM/MM on the fly frequencies using the VB (crosses), TM (squares), and GO (triangles) PPs 
with cut-off of 30 Ry, 70 Ry, and 140 Ry, respectively^ The solid line is a fit of experimental 
data according to Eq. ([3]) with lo^I2-kc = 3436 cm~^ [40], gq = 7.0911, ai = —5.7941 A , and 



a2 = 1.2711 



(b) Correlation between HB length N • • • O and fundamental N-H transition 



frequency (bullets: spectroscopic data for crystals contain ing intermolecular N-H • • • O HBs, solid 
line: fit according to Eq. ^ with Woo/2vrc = 3434 cm'^ \a% oq = -0.5345, oi = -0.8332 A~\ 
and 02 = 0.5044 A"^). 



the contribution from U:U dimers should be negligible 4l|]. A similar conclusion can be 
reached based on the work of Miller et al. who demonstrated an 1:1 stoichiometry for 
similar purine and pyrimidine base analogues in CDCI3 by monitoring the IR absorption 



42|. 



changes around 3210 cm~^ as a function of the uracil mole fraction 

Giving this evidence for the correct experimental assignment, we compare the present 
findings with recent gas phase quantum chemical studies of the thymine NH stretching vi- 
bration in the adenine: thymine base pair jsj. In harmonic approximation the DFT/B3LYP 
(6-31++G(d,p)) level of theory yields a frequency of 2981 cm~^ . Accounting for anhar- 
monicity within a three-dimensional model which includes all H-bonded stretching vibra- 
tions lowers this value to 2608 cm~^ . In order to rule out that this is a result of too soft 
DFT potential energy surfaces, MP2 corrections have been added to the diagonal anhar- 
monicity. However, this causes an increase of the thymine stretching frequency to 2688 
cm~^ only. Very recently, Yagi and coworkers have started to investigate this vibration us- 
ing their VMP2-(4) method |43j on the basis of a B3LYP/6-31G++(d,p) full-dimensional 



quartic force field potential. They obtained a preliminary value as low as 2488 cm ^ 44 1. 



At this level of theory the N ■ ■ ■ N distance is 2.88 A, which is smaller than the aver- 
age value obtained from the QM/MM trajectory. Therefore, the lower value of the NH 
stretching frequency in the gas phase is consistent with the stronger HB as compared to 
the solution phase. However, giving the empirical correlations discussed in the following 
the absolute value of the calculated frequencies appears to be considerably too small. 
To cope with this problem we decided to resort to an empirical mapping of HB distances 



to the stretching frequencies. In Ref. 2J] we have shown that this correlation plotted 
in Fig. [6^ gives an average frequency of 3227 cm^^ for the VB PP with a cut-off of 30 
Ry and using a 6.2 ps trajectory. Furthermore, the lineshape derived from the frequency 
fluctuations has been in good agreement with the experiment as well. This supports the 
assumption that despite the discrepancy in the absolute value of the stretching frequency 
the current QM/MM protocol provides a reasonable description of the fluctuations influ- 
encing the HB. We have extended this approach to the case of the N-H ■ ■ ■ O bond; the 
frequency-HB length correlation plot derived from available crystal data is shown in Fig. 
[6b. Notice that the data used included HBs involving amine groups, that is, the fact that 
the considered vibrations is actually of symmetric NH2 character is taken care of by the 
empirical relation. The resulting time-dependent frequencies for the two HBs are given 
in Fig. [71 
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FIG. 7: Time-dependent fundamental transition frequencies for the N-H • • • N (a) and the 
N-H • • • O (b) HB as obtained using the empirical correlations displayed in Fig. [6) The dot- 
ted lines are the averaged values {uJio) which are equal to 3132 cm^^ in (a) and 3275 cm^^ in 
(b). The dashed lines are the N • • • N and N • • • O distances scaled and shifted such as to reveal 
the correlation with the time-dependent frequencies. 
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FIG. 8: Frequency correlation function, Eq. ^ for the NH stretching vibrations of (a) the 
N-H • • • N and (b) the N-H • • • O HB. 

B. IR Lineshape Analysis 

The correlation function C (t) has been calculated for both HBs for the gap fluctuations 
shown in Fig. [3 The results given in Fig. [8] reveal a rapid initial decay followed by 
pronounced oscillations lasting up to the 2 ps which can be covered by the available 



trajectory (for a convergence study, see also Ref. [2J|)- In order to gain further insight 
into the distribution of modes coupled to the considered transitions, we have fitted C{t) 
according to the Brownian oscillator model, Eq. ([7]). The chosen frequency range of 0-1770 
cm~^ covers the low-frequency part of the calculated gas phase harmonic IR spectrum of 
the A:U pair. It is discretized into 20 modes whose frequencies and couphng strengths are 
used as fitting parameters. The distribution of coupling strengths S{uj) is shown in Fig. 
[HI First, we notice that both distributions are dominated by low-frequency vibrations. 
Second, in both cases we observe two smaller peaks which are around 200 and 500 cm~^ for 
N-H ■ ■ ■ N as well as around 350 and 550 cm~^ for N-H ■ ■ ■ O. It is tempting to assign these 
peaks to intermolecular vibrational modes. Indeed gas phase harmonic analysis of the 
A:U pair gives several intermolecular modes modulating the N ■ ■ ■ N and N ■ ■ ■ O distances 
in region of 150-400 cm~^ . Modes in the range of 400-600 cm~^ are mostly modulating 
the N ■ ■ ■ O distance consistent with Fig. [9l 

From the correlation functions one obtains the IR lineshape which is exemplarily shown 
for the N-H ■ ■ ■ N case in Fig. [TOl First we notice that the overall agreement with the 
experimental data is fairly good given the simple empirical model for assigning transition 
frequencies. In particular the asymmetry of the line is well reproduced. We also show the 



previous result obtained using the VB PP [24] which merely gives a slightly lower value 
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FIG. 9: Coupling strength between the NH vibration and the bath as function of the bath 
frequency, (a) results for N-H • • • N HB; (b) results for N-H • • • O HB. Points are fitted results 
and line is guide for eye. 

for the FHWM. 
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FIG. 10: IR lineshape of NH stretching compared with experimental data [g| (line for experi- 
mental result is drawn as guide for the eye only). Note that we have not included the shift of 
the transition frequency due to the imaginary part of the correlation function. It amounts to -3 
cm~^ for the TM and -16 cm~^ for the VB case. 



A detailed comparison of lineshape parameters obtained by using different PPs and 
cut-off parameters is presented in Tab. [H If we take the average transition frequencies and 
the corresponding energy gap between the two stretching vibrations as the measure for the 
performance of the method, we find that the TM PP with a cut-off of 70 Ry gives the best 
results. Here the deviations from experimental values are below 2 % for both frequencies. 
Of course, this result is not general in the sense that it is based on the proposed simulation 
protocol which includes the empirical mapping of distances to frequencies. Therefore 



TABLE I: Comparison of HB geometry and lineposition and lineshape for different PPs and 
cut-off parameters. Labels Ln...n/0) 2^nh---n/o ^n---n/o stand for bond length, fundamental 
transition frequency and the FHWM of the IR lineshape for NH • • • N/0 HB, respectively. Gap 
refers to the distance between the fundamental transition frequency of N-H • • • N and N-H • • • O 
HBs and L^^ .j,^ to the FHWM of the N-H • • • N HB based on the stretching-distance correlation 
from snapshot potentials. The TM (70 Ry) setup has been used in Fig. [TOl Note that the 
frequencies do not contain the shift due to the imaginary part of the correlation function which 



is smaller t 
from Ref. 



lan 5 cm in all cases. Further notice that the VB 30 Ry results are slightly different 



24| where a 6.2 ps trajectory had been used. 
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using the larger cut-off of 100 Ry does not necessarily give an improvement. Another 
point to mention is that Tab. [T] also shows that the correlation between average HB 
distance and average stretching frequency is not linear and in particular the fluctuations 
are sensitive to the used CPMD setup. The latter are also responsible for the linewidths 
which are compared in the right part of Tab. [B Also for this quantity the TM PP (70 
Ry) performs reasonably well for the case of the z/n_h - n transition. The FWHM of the 
i^N-H --o transition is, however, considerably underestimated in all cases. Finally, we give 
the FWHM obtained by using the gap correlation function derived from the snapshot 
potentials. They are remarkably close the values calculated via the empirical correlation. 
This leads us to conclude that the quantitative difference between the two approaches is 
indeed merely a constant shift of the transition frequencies (cf. Fig. [6]). 



IV. SUMMARY 



We have presented a QM/MM approach to the simulation of the intermolecular HB 
dynamics in solvated A:U pairs. On average the two HBs are almost linear and follow 
the empirical geometric relation of Limbach and coworkers between the HB lengths and 
the proton position. Further it was found that nonplanarity of the HBs is closely related 
to the compression due to bending motions of the C=0 and NH groups. The trajectory 
mean of the dihedral angle formed by HB donors and acceptors is predicted as -19°. Ex- 
ploring the empirical correlation between HB lengths and transition frequencies of the 
H-bonded protons we were able to simulate the IR spectrum in the region of the uracil 
NH and the adenine symmetric NH2 stretching fundamental transition. Attempts to use 
transition frequencies derived from snapshot potentials gave considerably red-shifted tran- 
sitions independent on the used pseudopotential and cut-off parameters. The difference 
between empirical and snapshot-derived frequency-distance correlation curves amounts to 
a shift of 450 cm~^ approximately independent on the HB distance. For the used empiri- 
cal mapping the TrouUier- Martins pseudopotential with a cut-off of 70 Ry gives the best 
agreement for the two fundamental transition. A multimode Brownian oscillator analysis 
of the gap correlation function yielded a spectrum of coupling strengths which is dom- 
inated by low-frequency vibrations, but contains contributions presumably from modes 
which modulate the two HBs in the range between 200-700 cm~^ . While the width of 
the i^NH - N transition could be fairly well reproduced, the width of the i^NH - o transition is 
underestimated. However, one should be aware that our simulation does neither include 
contributions from other high-frequency modes nor combination bands which could give 
rise to an additional broadening of the experimental spectra. 
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